QTW Case Study 4 - ARIMA Modelling using FluNet data from the World Health Organization

Jayson Barker, Brandon Croom, Shane Weinstock

Business Understanding

The World Health Organization (WHO) has created a publically accessible repository of data specific to Flu prevalance called FluNet. The data providers are all collaborating laboratories / national focal points from different countries, areas or territories networked with the WHO. FluNet is a global web-based tool for influenza virological surveillance first launched in 1997. The data at country level are publically available and updated weekly. The data are provided by National Influenza Centres (NICs) of the Global Influenza Surveillance and Response System (GISRS) and other national influenza reference laboratories collaborating actively with GISRS, or are uploaded from WHO regional databases. Using this data source, our team identified an appropriate data source that contained influenza type A and B prevalance by week and year to perform ARIMA modelling on. The weekly data is well suited for ARIMA testing as it provides us weekly set of flu data that is well suited for ARIMA.

An autoregressive integrated moving average (ARIMA) model is one of the many models used in time series analysis. ARIMA models are applied in cases where data show evidence of non-stationarity. Non-stationarity typically occurs when a the data points of the process being modeled are dependent upon time. An ARMINA model is essentially a combination of AR and MA models. The AR part of ARIMA indicates that the variable of interest is regressed on its own prior values. The AR model order is expressed by the value p. The MA part indicates that the regression error is a linear combination of error terms whos values occurred at various times in the past. The MA model order is expressed by the value q. The I indicates that the data values have been differenced and indicates the number of differences. The I order of the model is indicated by the value d. Typically ARIMA models are denoted as ARIMA(p,d,q), where the p value represents the order of the AR model portion, d represents the order of the I portion of the model, and the q value represents the order of the MA model portion.

Another model that is used in against data that is seasonal is the seasonal auto regressive integrated moving average with exogenous factors model, also known as the SARIMAX model. SARIMAX modeling is similar to ARIMA modeling and maintains the same types of p, d, and q values for AR, I and MA respectively. SARIMAX also introduces an s component, which indicates the seasonality of the data. The difference between ARIMA and SARIMAX is the seasonality and exogenous factors.

The parameters we used to extract this data include all influenza type A & B records from January 2016 through September 2020. This provides us a robust data set that includes influenza type A & B positive case counts as well as negative case counts. The overall goal of this report is to show that influenza data can be leveraged within an ARIMA model.

# import libraries
import pandas as pd
import pandas_profiling
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.tsa.stattools import adfuller
from numpy import log
%matplotlib inline
from matplotlib import pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pmd
from pmdarima.arima.utils import ndiffs
#from pmdarima.arima import ADFTest

# Define helper functions
def arimamodel(timeseriesarray, test):
    autoarima_model = pmd.auto_arima(timeseriesarray, 
                              start_p=0, 
                              start_q=0,
                              test=test,
                              trace=True)
    return autoarima_model

# Define settings values for plotting
plt.rcParams.update({'figure.figsize':(15,10), 'figure.dpi':120})

Exploratory Data Analysis

In this section we will discuss the steps and results found as part of the exploratory data analysis (EDA). During EDA we investigate the data looking for data anomalies, missing data and any other data work that needs to occur prior to model building. We first start by loading the data that was extracted from the FluNet website. The data is focused on the US only and covers the 2016 - 2020 timespan

# load data into dataframe
df = pd.read_csv('FluNetInteractiveReport.csv')
#df.columns.tolist()

Leveraging the Pandas Profiling tool provides a wealth of information on the data set and helps to speed up the EDA. Summary observations from pandas profiling include:

- There are 22 variables

- There are 246 observations

- Only 5.2% of the data set has missing values

Digging into the specific fields that will be required for our analysis. There are key fields that we must pay special attention to: ALL_INF, SDATE and EDATE.

- ALL_INF is our column of interest as it includes the number of positive cases identified during that specific week of the year.

- SDATE is the start date of the week the data was reported

- EDATE is the end date of the week the data was reported

For these three fields pandas profiling indicates the follwing information:

- SDATE has no missing values, is 100% distinct and is a categorical variable

- EDATE has no missing values, is 100% distinct and is a categorical variable

- ALL_INF has no missing values, is 93.5% distinct and is a numerical variable

All of this information allows us to understand the following about the data:

- There are no duplicate weeks that need to be addressed

- There is no missing data in the fields required for analysis that need to be addressed

- SDATE and EDATE will need to be converted to datatime data types to perform time series analysis

# Profile
profile = df.profile_report(title = "Flu Dataset")
profile.to_notebook_iframe()

# Number of specimes columns: Received/Coll, Processed
# Number of influenza A viruses detected by subtype: A(H1), A (H1N1) pdm09, A (H3), A (H5), A (not subtyped), A (total)
# Number of influenza B viruses detected by subtype: B (Yamagata), B (victoria), B (lineage not determined), B (total)
# Total number of influenza positive
# Total number of influenza negative
# ILI activity

Summarize dataset: 100%|██████████| 36/36 [03:32<00:00,  5.90s/it, Completed]
Generate report structure: 100%|██████████| 1/1 [00:43<00:00, 43.49s/it]
Render HTML: 100%|██████████| 1/1 [00:22<00:00, 22.44s/it]

Below is a plot of the original data set to show the basic data trends over time.

# Convert the EDATE and SDATE values to datetime values
df['EDATE'] = pd.to_datetime(df['EDATE'])
df['SDATE'] = pd.to_datetime(df['SDATE'])

#build a second data frame for ADF testing and autoarima testing
df2 = df[['EDATE','ALL_INF']]
df2.set_index('EDATE',inplace=True)

#plot the data set
df2.plot()

<AxesSubplot:xlabel='EDATE'>

ARIMA Analysis

In this section we will begin the ARIMA analysis. In looking at the plot of the data shown above, we see what looks to be seasonal values. This would lead us to believe the data may not be stationary. Stationarity can be tested in different ways and as such we'll look at three differnt statiionarity tests to get a determination on stationarity. We'll look at three different tests:

- Augmented Dickey-Fuller (ADF) Test

- Phillips-Peron (PP) Test

- Kwiatkowski=Phillips=Schmidt-Shin (KPSS) Test

adf_test = pmd.arima.ADFTest(alpha=0.05)
adf_test.should_diff(df2)

(0.01, False)

n_adf_adf = pmd.arima.utils.ndiffs(df2, test='adf')
n_adf_pp = pmd.arima.utils.ndiffs(df2, test='pp')
n_adf_kpss = pmd.arima.utils.ndiffs(df2, test='kpss')

print("ADF Test Results: ", n_adf_adf)
print("PP Test Results: ", n_adf_pp)
print("KPSS Test Results: ", n_adf_kpss)

ADF Test Results:  0
PP Test Results:  1
KPSS Test Results:  0

The results from the stationarity tests seem to indicate the data is non-stationary. The ADF and the KPSS test is detecting that the data is stationary, while the PP method is indicating that the data is non-stationary. Given the domain knowledge available on flu data, it would seem that treating the data as non-stationary would be a good first test.

In the plots below, we perform our first run of auto-correlation to evaluate lags and what parameters may make the most sense for our ARIMA model. In the original series, we see peaks left plot that seem to indicate a seasonal type pattern. The autocorrelation plot shows a cyclic like pattern that typically is representative of seasonal values. In order to remove the seasonality from the data, differencing is needed. Differencing is a method of transforming a time series analysis to remove dependence on time. The second plot below shows the first difference being taken on the data. In the left plot we see some of the seasonality, but the data is starting to move more toward a less time dependent mean. Differencing a second time appears to have the best alignment to the mean, which sets our lag at 3 - right above the cutoff per the autocorrelation visual.

# Original Series
fig, axes = plt.subplots(3, 2, sharex=True)
axes[0, 0].plot(df['ALL_INF']); axes[0, 0].set_title('Original Series')
plot_acf(df['ALL_INF'], ax=axes[0, 1])

# 1st Differencing
axes[1, 0].plot(df['ALL_INF'].diff()); axes[1, 0].set_title('1st Order Differencing')
plot_acf(df['ALL_INF'].diff().dropna(), ax=axes[1, 1])

# 2nd Differencing
axes[2, 0].plot(df['ALL_INF'].diff().diff()); axes[2, 0].set_title('2nd Order Differencing')
plot_acf(df['ALL_INF'].diff().diff().dropna(), ax=axes[2, 1])

# Stationary at 2nd difference

plt.show()

Here we perform our auto-correlation testing for additional ARIMA parameters.

# PACF plot of 1st differenced series
plt.rcParams.update({'figure.figsize':(9,3), 'figure.dpi':120})

fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].plot(df['ALL_INF'].diff().diff()); axes[0].set_title('2nd Differencing')
axes[1].set(ylim=(0,5))
plot_pacf(df['ALL_INF'].diff().dropna(), ax=axes[1])

plt.show()

In this section, we evaluate auto-correlation plot to determine the q value.

# MA term

fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].plot(df['ALL_INF'].diff()); axes[0].set_title('1st Differencing')
axes[1].set(ylim=(0,1.2))
plot_acf(df['ALL_INF'].diff().dropna(), ax=axes[1])

# q = 3

plt.show()

Using the values defined above for P, D, and Q, we run our ARIMA model with values, 1, 1, and 3. This identifies that the constant is significant, and the moving average is not closed to 0. We can also see that the pvalue of AR and MA terms are highly significant at < 0.05, which indicates this model is sound.

# 1,1,3 ARIMA Model
# p d q
model = ARIMA(df['ALL_INF'], order=(1,1,3))
model_fit = model.fit(disp=0)
print(model_fit.summary())

# note constant is significant
# ma is not close to zero
# pvalue of AR and MA terms are highly significant < 0.05

                             ARIMA Model Results                              
==============================================================================
Dep. Variable:              D.ALL_INF   No. Observations:                  245
Model:                 ARIMA(1, 1, 3)   Log Likelihood               -2067.006
Method:                       css-mle   S.D. of innovations           1114.296
Date:                Mon, 19 Oct 2020   AIC                           4146.012
Time:                        21:48:25   BIC                           4167.020
Sample:                             1   HQIC                          4154.472
                                                                              
===================================================================================
                      coef    std err          z      P>|z|      [0.025      0.975]
-----------------------------------------------------------------------------------
const              -3.8501    172.387     -0.022      0.982    -341.722     334.022
ar.L1.D.ALL_INF    -0.4297      0.167     -2.566      0.010      -0.758      -0.101
ma.L1.D.ALL_INF     1.2121      0.161      7.524      0.000       0.896       1.528
ma.L2.D.ALL_INF     0.8519      0.141      6.040      0.000       0.575       1.128
ma.L3.D.ALL_INF     0.4105      0.071      5.779      0.000       0.271       0.550
                                    Roots                                    
=============================================================================
                  Real          Imaginary           Modulus         Frequency
-----------------------------------------------------------------------------
AR.1           -2.3272           +0.0000j            2.3272            0.5000
MA.1           -1.2639           -0.0000j            1.2639           -0.5000
MA.2           -0.4056           -1.3277j            1.3882           -0.2972
MA.3           -0.4056           +1.3277j            1.3882            0.2972
-----------------------------------------------------------------------------

The residuals plot below shows the distribution of the mean about 0 - and also does not indicate a shift in the mean (i.e., a constant mean is present) which gives us additional confidence in our ARIMA model.

# Plot residual errors
residuals = pd.DataFrame(model_fit.resid)
fig, ax = plt.subplots(1,2)
residuals.plot(title="Residuals", ax=ax[0])
residuals.plot(kind='kde', title='Density', ax=ax[1])
plt.show()

# we see a constant mean and no identifiable patterns

Fitting the actual data onto our predicted ARIMA models demonstrates a good fit. This shows that the ARIMA(2,0,0) model is a good approximation of the model.

# Actual vs Fitted
model_fit.plot_predict(dynamic=False)
plt.show()

# dynamic=False causes the lags to be used for prediction

Another option to determine the p,d,q values of the ARIMA model is to leverage auto ARIMA functionality. With auto ARIMA the intent is to execute multiple models with multiple values of p,d,q to determine the best model configuration for the data. The output of the auto ARIMA execution is shown below. The auto ARIMA function allows for testing against the same stationarity testing models to build out an ARIMA model. For purposes of determining the best model using auto ARIMA we'll leverage all three tests

arima_model_adf = arimamodel(df2,"adf")
arima_model_adf.summary()

Performing stepwise search to minimize aic
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=5030.440, Time=0.07 sec
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=4320.770, Time=0.12 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=4724.586, Time=1.00 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=5119.955, Time=0.02 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=4156.330, Time=0.14 sec
 ARIMA(3,0,0)(0,0,0)[0] intercept   : AIC=4158.075, Time=0.22 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=4158.052, Time=0.26 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=4204.781, Time=0.67 sec
 ARIMA(3,0,1)(0,0,0)[0] intercept   : AIC=4159.284, Time=0.48 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=4160.351, Time=0.11 sec

Best model:  ARIMA(2,0,0)(0,0,0)[0] intercept
Total fit time: 3.114 seconds

SARIMAX Results

Dep. Variable:	y	No. Observations:	246
Model:	SARIMAX(2, 0, 0)	Log Likelihood	-2074.165
Date:	Mon, 19 Oct 2020	AIC	4156.330
Time:	21:48:34	BIC	4170.351
Sample:	0	HQIC	4161.976
	- 246
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
intercept	208.6941	169.694	1.230	0.219	-123.901	541.289
ar.L1	1.6506	0.027	61.886	0.000	1.598	1.703
ar.L2	-0.6988	0.026	-27.286	0.000	-0.749	-0.649
sigma2	1.221e+06	5.25e+04	23.230	0.000	1.12e+06	1.32e+06

Ljung-Box (Q):	52.82	Jarque-Bera (JB):	582.29
Prob(Q):	0.08	Prob(JB):	0.00
Heteroskedasticity (H):	2.13	Skew:	0.13
Prob(H) (two-sided):	0.00	Kurtosis:	10.53

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

arima_model_adf = arimamodel(df2,"pp")
arima_model_adf.summary()

Performing stepwise search to minimize aic
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=4302.494, Time=0.04 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=4152.832, Time=0.07 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=4188.281, Time=0.34 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=4300.497, Time=0.02 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=4152.191, Time=0.11 sec
 ARIMA(3,1,0)(0,0,0)[0] intercept   : AIC=4154.137, Time=0.14 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=4154.187, Time=0.27 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=4152.429, Time=0.12 sec
 ARIMA(3,1,1)(0,0,0)[0] intercept   : AIC=4156.070, Time=0.45 sec
 ARIMA(2,1,0)(0,0,0)[0]             : AIC=4150.193, Time=0.06 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=4150.834, Time=0.05 sec
 ARIMA(3,1,0)(0,0,0)[0]             : AIC=4152.138, Time=0.10 sec
 ARIMA(2,1,1)(0,0,0)[0]             : AIC=4152.188, Time=0.19 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=4150.430, Time=0.08 sec
 ARIMA(3,1,1)(0,0,0)[0]             : AIC=4154.072, Time=0.21 sec

Best model:  ARIMA(2,1,0)(0,0,0)[0]          
Total fit time: 2.328 seconds

SARIMAX Results

Dep. Variable:	y	No. Observations:	246
Model:	SARIMAX(2, 1, 0)	Log Likelihood	-2072.097
Date:	Mon, 19 Oct 2020	AIC	4150.193
Time:	21:48:36	BIC	4160.697
Sample:	0	HQIC	4154.423
	- 246
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
ar.L1	0.7490	0.043	17.486	0.000	0.665	0.833
ar.L2	-0.1034	0.050	-2.075	0.038	-0.201	-0.006
sigma2	1.308e+06	5.53e+04	23.660	0.000	1.2e+06	1.42e+06

Ljung-Box (Q):	44.74	Jarque-Bera (JB):	649.67
Prob(Q):	0.28	Prob(JB):	0.00
Heteroskedasticity (H):	2.08	Skew:	-0.95
Prob(H) (two-sided):	0.00	Kurtosis:	10.75

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

arima_model_adf = arimamodel(df2,"kpss")
arima_model_adf.summary()

Performing stepwise search to minimize aic
 ARIMA(0,0,0)(0,0,0)[0] intercept   : AIC=5030.440, Time=0.03 sec
 ARIMA(1,0,0)(0,0,0)[0] intercept   : AIC=4320.770, Time=0.09 sec
 ARIMA(0,0,1)(0,0,0)[0] intercept   : AIC=4724.586, Time=0.59 sec
 ARIMA(0,0,0)(0,0,0)[0]             : AIC=5119.955, Time=0.02 sec
 ARIMA(2,0,0)(0,0,0)[0] intercept   : AIC=4156.330, Time=0.11 sec
 ARIMA(3,0,0)(0,0,0)[0] intercept   : AIC=4158.075, Time=0.20 sec
 ARIMA(2,0,1)(0,0,0)[0] intercept   : AIC=4158.052, Time=0.20 sec
 ARIMA(1,0,1)(0,0,0)[0] intercept   : AIC=4204.781, Time=0.65 sec
 ARIMA(3,0,1)(0,0,0)[0] intercept   : AIC=4159.284, Time=0.47 sec
 ARIMA(2,0,0)(0,0,0)[0]             : AIC=4160.351, Time=0.13 sec

Best model:  ARIMA(2,0,0)(0,0,0)[0] intercept
Total fit time: 2.505 seconds

SARIMAX Results

Dep. Variable:	y	No. Observations:	246
Model:	SARIMAX(2, 0, 0)	Log Likelihood	-2074.165
Date:	Mon, 19 Oct 2020	AIC	4156.330
Time:	21:48:39	BIC	4170.351
Sample:	0	HQIC	4161.976
	- 246
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
intercept	208.6941	169.694	1.230	0.219	-123.901	541.289
ar.L1	1.6506	0.027	61.886	0.000	1.598	1.703
ar.L2	-0.6988	0.026	-27.286	0.000	-0.749	-0.649
sigma2	1.221e+06	5.25e+04	23.230	0.000	1.12e+06	1.32e+06

Ljung-Box (Q):	52.82	Jarque-Bera (JB):	582.29
Prob(Q):	0.08	Prob(JB):	0.00
Heteroskedasticity (H):	2.13	Skew:	0.13
Prob(H) (two-sided):	0.00	Kurtosis:	10.53

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

The output of the auto ARIMA shows that an ARIMA(2,0,0) is deemed the best model when using the ADF and KPSS stationarity test. This is based off of the AIC value of 4156.330 being the lowest of the models executed. The auto ARIMA output using the PP stationarity test was an ARIMA(2,1,0). Logically this model selection makes sense as we know there is seasonality in the data. In the output it is also seen that the SARIMAX(2,1,0) model looks like a viable model and has the same AIC as the ARIMA model. This helps to reinforce the similarities between ARIMA and SARIMAX models. Between the three models the ARIMA(2,1,0) also has the lowest AIC at 4150.193. Not a huge difference from the other two but lower none the less.

In order to test the forecasting of the ARIMA we will create a train/test split of the data and reexecute the model. The predictions from the model will become the forecast. This approach will allow us to see how close the model forecasts. Below we'll build out our train/test split and plot the training and test splits to visualize which part of the data the forecast will go against.

train = df2[:216]
test = df2[-30:]
plt.plot(train, label="Training")
plt.plot(test, label='Test')
plt.legend(loc='upper left')

<matplotlib.legend.Legend at 0x148e0378a30>

Next we'll execute the auto ARIMA against the training data set to determine which model is best suited for this analysis. As we can see from the model output below, as with the full data set the auto ARIMA selected the ARIMA(2,0,0).

arima_model = arimamodel(train,"pp")
arima_model.summary()

Performing stepwise search to minimize aic
 ARIMA(0,1,0)(0,0,0)[0] intercept   : AIC=3764.918, Time=0.03 sec
 ARIMA(1,1,0)(0,0,0)[0] intercept   : AIC=3647.291, Time=0.07 sec
 ARIMA(0,1,1)(0,0,0)[0] intercept   : AIC=3676.863, Time=0.34 sec
 ARIMA(0,1,0)(0,0,0)[0]             : AIC=3763.739, Time=0.03 sec
 ARIMA(2,1,0)(0,0,0)[0] intercept   : AIC=3647.768, Time=0.11 sec
 ARIMA(1,1,1)(0,0,0)[0] intercept   : AIC=3648.074, Time=0.11 sec
 ARIMA(2,1,1)(0,0,0)[0] intercept   : AIC=3649.569, Time=0.42 sec
 ARIMA(1,1,0)(0,0,0)[0]             : AIC=3645.372, Time=0.04 sec
 ARIMA(2,1,0)(0,0,0)[0]             : AIC=3645.883, Time=0.06 sec
 ARIMA(1,1,1)(0,0,0)[0]             : AIC=3646.176, Time=0.07 sec
 ARIMA(0,1,1)(0,0,0)[0]             : AIC=3675.276, Time=0.15 sec
 ARIMA(2,1,1)(0,0,0)[0]             : AIC=inf, Time=0.53 sec

Best model:  ARIMA(1,1,0)(0,0,0)[0]          
Total fit time: 1.987 seconds

SARIMAX Results

Dep. Variable:	y	No. Observations:	216
Model:	SARIMAX(1, 1, 0)	Log Likelihood	-1820.686
Date:	Mon, 19 Oct 2020	AIC	3645.372
Time:	21:48:42	BIC	3652.114
Sample:	0	HQIC	3648.096
	- 216
Covariance Type:	opg

	coef	std err	z	P>|z|	[0.025	0.975]
ar.L1	0.6578	0.027	23.974	0.000	0.604	0.712
sigma2	1.332e+06	6.07e+04	21.949	0.000	1.21e+06	1.45e+06

Ljung-Box (Q):	40.19	Jarque-Bera (JB):	521.77
Prob(Q):	0.46	Prob(JB):	0.00
Heteroskedasticity (H):	2.05	Skew:	-0.82
Prob(H) (two-sided):	0.00	Kurtosis:	10.45

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

prediction = pd.DataFrame(arima_model.predict(n_periods=30),index=test.index)
prediction.columns = ['ALL_INF']

Utilizing the prediction capability of the model we can test to see how well the model fits the test portion of the data. The plot below shows the output of those results. From the plot we see the prediction mimics the test set fairly well. It follows the downward slope but does begin to trend back up earlier than the test data.

plt.plot(train,label="Training")
plt.plot(test,label="Test")
plt.plot(prediction,label="Prediction")
plt.legend(loc='upper left')
plt.show()

Conclusion

In conclusion, leveraging the flu data for the US from 2016-2020 we have been able to validate that time series modeling can be performed on the data and achieve reasonable results. Over the course of this report we have shown the following items:

- Check for stationarity within a data set to determine how to approach modeling

- Determining manually the best p,d,q values for an ARIMA model

- Forecasting of the manual model

- Determining through auto ARIMA the best p,d,q values for an ARIMA model

- Forecasting of the auto ARIMA model

Over the course of the report both for the manual and auto ARIMA there were differences in the ARIMA models developed. The manual ARIMA model was an ARIMA(1,1,3) while the best auto ARIMA model resulted in an ARIMA(2,1,0). The differences in these models further manifested in the forecasts that were generated. The ARIMA(1,1,3) model did a much better job of forecasting than the ARIMA(2,0,0). The ARIMA(1,1,3) forecast followed very closely the actual data values with its forecast. The ARIMA(2,1,0) initially followed the forecast well but then deviated quickly from the test values of the model. With respect to AIC values the manual ARIMA(1,1,3) model had an AIC of 4146.01, while the auto ARIMA(2,1,0) model had an AIC of 4150.193. The manual model ARIMA still outperformed the auto ARIMA model in AIC as well. All in all, the models generated had the ability to forecast seasonal flu trends and were able to perform pretty well